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Abstract. Numerical methods for solving the ideal magnctohydrodynamic (MHD) equations in more than 
one space dimension must confront the challenge of controlling errors in the discrete divergence of the magnetic 
field. One approach that has been shown successful in stabilizing MHD calculations are constrained transport 
(CT) schemes. CT schemes can be viewed as predictor-corrector methods for updating the magnetic field, 
where a magnetic field value is first predicted by a method that does not exactly preserve the divergence-free 
condition on the magnetic field, followed by a correction step that aims to control these divergence errors. In 
I * Hclzcl ct al. [J. Comp. Phys. 227, 9527 (2011)] the authors presented an unstaggcrcd constrained transport 

method for the MHD equations on three-dimensional Cartesian grids. In this approach an evolution equation 

(^*\ for the magnetic potential is solved during each time step and a divergence-free update of the magnetic field 

is computed by taking the curl of the magnetic potential. The evolution equation for the vector potential is 

l/~) only weakly hyperbolic, which requires special numerical treatment. A key step in this method is the use of 

^_^ dimensional splitting in order to overcome these difficulties. 

In this work we generalize the method of [J. Comp. Phys. 227, 9527 (2011)] in three important ways: (1) 

j ' we remove the need for operator splitting by switching to an appropriate method of lines discretization and 

^L^ coupling this with a non-conservative finite volume method for the magnetic vector potential equation, (2) we 

^~^ increase the spatial and temporal order of accuracy of the entire method to third order, and (3) we develop 

the method so that it is applicable on both Cartesian and logically rectangular mapped grids. The method 

of lines approach that is used in this work is based on a third-order accurate finite volume discretization in 

space coupled to a third order strong stability preserving Rungc-Kutta time stepping method. The evolution 

equation for the magnetic vector potential is solved using a non-conservative finite volume method based on the 

approach of Castro et al. [Math. Cornput. 79, 1427 (2010)]. The curl of the magnetic potential is computed 

via a third-order accurate discrete operator that is derived from appropriate application of the divergence 

theorem and subsequent numerical quadrature on element faces. Special artificial resistivity limiters are used 

£f) to control unphysical oscillations in the magnetic potential and field components across shocks. Several test 

^. computations are shown that confirm third order accuracy for smooth test problems and high-resolution for 

f*"*) test problems with shock waves. 
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1. Introduction. It is well known that numerical methods for the multi-dimensional 
magnctohydrodynamic (MHD) equations must in some way respect the divergence-free con- 
dition on the magnetic field; failure to control divergence errors have been shown to lead to 
nonlinear numerical instabilities (e.g., see Toth [42] for a discussion). Brackbill and Barnes 
[9] were one of the first to document this problem; since their paper a variety of approaches 
have been proposed to control the divergence of the magnetic field. The four main approaches 
that have been developed are 8- wave [36, 37], projection [3, 42], divergence-cleaning [18], and 
constrained transport methods [2, 4, 17, 20, 21, 26, 30, 31, 39, 40, 41, 42, 45]. We will not 
endeavor to provide a detailed description of these four numerical approaches, since such a 
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description already exists in Toth [42] and Hclzcl et al. [26]. Instead, because in this work we 
will develop a constrained transport method, we very briefly describe the CT methodology 
for MHD in the remainder of this section. 

The basic framework of the constrained transport methodology was developed by Evans 
and Hawley [20]. The premise was to modify the celebrated Yee scheme [48] from vacuum 
electromagnetism to the ideal MHD equations; this framework included staggered electric 
and magnetic fields, just as in the Yee scheme. The resulting method can be viewed as 
predictor-corrector approach for updating the magnetic field, where a magnetic field value is 
first predicted by a method that does not exactly preserve the divergence-free condition on 
the magnetic field, followed by a correction step that aims to control these divergence errors. 
The basic steps in the Evans and Hawley CT approach can be outlined as follows: 

Step 0. Start with MHD cell average values at the current time: (p n , pu n , £ n , B"). 

Step 1. Take a time-step using some finite volume method, which produces the cell average 

values: (p n+1 , pu n+1 , £ n+1 , B*). B* is the predicted value of the magnetic field. 
Step 2. Using the ideal Ohm's law relationship, E = B x u, and some space and time 

interpolation scheme for B and u, reconstruct a space and time staggered electric 

field value: E™+5. 
Step 3. Produce the corrected magnetic field value, B™ +1 , from the Yee scheme [48] step for 

the magnetic field: 

B" +1 =B"-AiVx E n+ 5, 

where the curl operation is understood to be a discrete curl. 

An alternative, but completely equivalent, interpretation of this scheme replaces Step 3 above 
with the following: 

Step 3a. Produce the magnetic potential value, A™ +1 , from the induction equation written 
in potential form: 

A «+i = A™- AtE"+5. 

Step 3b. Produce the corrected magnetic field value, B™ +1 , by computing the discrete curl 
of the magnetic potential: 

B n+1 = y x A «+l_ 

By altering how Step 2 (i.e., the reconstruction of the electric field) is done in the Evans 
and Hawley framework, a variety of CT schemes were developed, including Balsara and Spicer 
[4], Dai and Woodward [17], Ryu et al. [40], and Londrillo and Del Zanna [30, 31]. In his 
review article, Toth [42] showed that constrained transport methods could also be formulated 
without the use of staggered electric and magnetic fields. Several of these unstaggered variants 
have been formulated, including the schemes of Fey and Torrilhon [21], Helzel et al. [26], and 
Rossmanith [39]. 

The particular focus of this work is to consider modifications of the constrained transport 
approaches of Rossmanith [39] and Helzel et al. [26]. Rossmanith [39] developed an unstag- 
gered constrained transport method for the two-dimensional MHD equations on Cartesian 
grids. In this approach LeVeque's wave propagation method [28, 29] is used to update the 
conserved variables of the MHD equations. Step 3a of the CT framework is modified by 
directly solving a scalar transport equation for the out-of-plane magnetic potential: 

A 3 t + u 1 A 3 x + u 2 A 3 y = 0. (1.1) 



The magnetic field components are computed via the curl operation on the magnetic potential: 



B 1 = A 3 and B 2 = -A 3 



(1.2) 



Transport equation (1.1) for the potential is solved using a modified version of LeVeque's wave 
propagation algorithm. In particular, the flux limiters of the wave propagation method are 
modified in the numerical solution of (1.1) in order to avoid not only unphysical oscillations 
in the magnetic potential but also in derivatives of the potential (i.e., the magnetic field). 

Hclzel et al. [26] extended this unstaggered constrained transport method to the three- 
dimensional MHD equations. In this case one has to consider the following transport equation 
for the vector potential: 
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(1.3) 



The magnetic field components are computed via the curl operation on the magnetic potential: 



B l =A 3 y -A 2 z , 



B- 



— Az , x , 



and B' 



— -™- tX A y . 



(1.4) 



In contrast to the 2D case, the 3D transport equation for the magnetic potential needs an ad- 
ditional constraint (i.e., the gauge condition) in order to get a closed set of evolution equations 
for the vector potential A and the scalar potential ip. Hclzel et al. [26] chose the so-called 
Weyl gauge (ip = 0) , which results in a transport equation for the vector potential that is only 
weakly hyperbolic. This means that standard methods for hyperbolic problems that are based 
on an eigenvector decomposition of the jump of the advected quantities at grid cell interfaces 
cannot be used. In order to overcome this difficulty, Hclzel et al. [26] developed a second 
order accurate dimensional split method for the magnetic vector potential equation that is 
well defined also in the weakly hyperbolic case. Many test computations confirmed that this 
approach leads to a robust and accurate approximation of the MHD equations on Cartesian 
grids. For smooth solutions of the MHD equations, the method from [26] was shown to be 
second order accurate. 

The focus of this work is to generalize the approach [26] in three important ways: 

1. We remove the need for operator splitting by switching to an appropriate method of 
lines discretization and coupling this with a non-conservative finite volume method 
for the magnetic vector potential equation. 

2. We increase the spatial and temporal order of accuracy of the entire method to third 
order using a third-order SSP-RK time-stepping method and a third-order finite vol- 
ume spatial discretization. 

3. We develop the method so that it is applicable on both Cartesian and logically rect- 
angular mapped grids. 

The nonconscrvativc method that we use on the magnetic vector potential equation is inspired 
by previous work on the approximation of non-conservative hyperbolic systems by Cancstrclli 
et al. [12], Dumbser et al. [19], and Ketcheson et al. [27]. By using a method of lines approach 
with a stong stability preserving Runge-Kutta (SSP-RK) method in time and third order ac- 
curate reconstruction in space, we are able to construct a third order accurate unstaggered 
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constrained transport method. In this approach, the correction of the magnetic field is per- 
formed during each stage of the Runge-Kutta method. For the spatial discretization we use a 
multidimensional (limited) polynomial reconstruction of the cell average quantities. 

After a brief review of the ideal MHD equations in Section 2, we introduce the proposed 
constrained transport scheme in a series of four sections. In Section 3 we outline the time 
discretization of the proposed method and show how to embed the constrained transport steps 
into the various stages of the third-order SSP Runge-Kutta scheme. In Section 4 we describe 
the third order accurate spatial discretization for the MHD system. In Section 5 we show how 
to modify this approach to obtain a high-order non-conservative spatial discretization of the 
magnetic vector potential equation. In this section we also develop special artificial resistivity 
limiters that are used to control unphysical oscillations in both the magnetic potential and 
magnetic field components. In Section 6 we show how to compute the curl to high-order of 
the discrete magnetic vector potential. Several numerical examples on Cartesian and mapped 
grids are shown in Section 7; these examples are used to both verify the third-order accuracy 
in space and time, as well as the shock-capturing ability of the proposed scheme. Conclusions 
can be found in Section 8. 

2. The ideal MHD equations. The ideal MHD equations are a first order hyperbolic 
system of conservation laws that can be written in the form 
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(2.1) 



(2.2) 



where p, pu and £ are the total mass, momentum and energy densities, and B is the magnetic 
field. The termal pressure, p, is related to the conserved quantities through the ideal gas law 



/>=(--l)(£-^l|B|| 2 -^||u|| 2 



(2.3) 



where 7 = 5/3 is the ideal gas constant. Here || • || denotes the Euclidean vector norm. A 
complete derivation and discussion of MHD system (2.1)-(2.2) can be found in several standard 
plasma physics textbooks (e.g., [14, 23, 33]). 

2.1. Hyperbolicity of the MHD system. We first note that system (2.1), along 
with the equation of state (2.3), provides a full set of equations for the time evolution of all 
eight state variables: (p, pu, £ , B). These evolution equations form a hyperbolic system. In 
particular, the eigenvalues of the flux Jacobian in some arbitrary direction n (||n|| = 1) can 
be written as follows: 



A 1 ' 8 = u • n =F Cf : fast magnetosonic waves, 
A 2 ' 7 = u • n q= c a : Alfven waves, 
A 3 ' 6 = u • n =F c s : slow magnetosonic waves, 
A = u • n : entropy wave, 



A 5 = u • n 



divergence wave, 

1 



(2.4) 
(2.5) 
(2.6) 
(2.7) 
(2.8) 
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The eigenvalues are well-ordered in the sense that 

A 1 < A 2 < A 3 < A 4 < A 5 < A 6 < A 7 < A 8 . 



(2.9) 
(2.10) 

(2.11) 
(2.12) 

(2.13) 



The fast and slow magnetosonic waves are genuinely nonlinear, while the remaining waves are 
linearly degenerate. Note that the so-called divergence-wave has been made to travel at the 
speed un via the Godunov-Powcll formulation [22, 36, 37], thus restoring Galilean invariance. 

2.2. The role of V • B = in numerical discretizations. If the initial magnetic field 
is divergence- free, then under the evolution described by equation (2.1), the divergence- free 
condition is satisfied for all time. This makes (2.2) an involution [16]. Unfortunately, although 
(2.2) is an involution under the evolution of the exact MHD system, it is generally not an 
involution of the discretized MHD equations. Furthermore, it has been well-documented in 
the literature that the failure of numerical methods to control divergence errors can lead to 
nonlinear numerical instabilities [42]. 

Explanations for why the divergence errors lead to numerical instabilities have been con- 
sidered from a few different points-of-view. We briefly review the two main explanations 
below. 

1. Brackbill and Barnes [8] showed that if V • B ^ 0, then the magnetic force, 

1 



BB 
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in the direction of the magnetic field, will not in general vanish: 

F B = ||B|| 2 (V-B)^O. 



(2.14) 



(2.15) 



If this spurious forcing becomes too large, it can lead to numerical instabilities [8, 39, 
42]. 
2. Barth [7] gave an explanation based on the well-known result of Godunov [22] that 
the MHD entropy density, 



U(q) = -plog(pp 7 ) 



(2.16) 



produces a set of entropy variables, v = U. q , that do not immediately symmetrize 
the ideal MHD equations. Instead, a symmetric hyperbolic form of ideal MHD can 
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only be obtained if an additional term that is proportional to the divergence of the 
magnetic field is included in the MHD equations: 



g,t + V-F(g) 



ideal MHD 



additional term 



o. 



where x{l) = (7 — 1) 



pu-B 



(2.17) 



By looking at how the entropy behaves on the discrete level, Barth [7] was able to 
prove that certain discontinuous Galerkin discretizations of the ideal MHD equtions 
could be made to be entropy stable (see Tadmor [43]) if the discrete magnetic field 
was made globally divergence-free. The implication of this result is that schemes that 
do not control errors in the divergence of the magnetic field run the risk of becoming 
entropy unstable. 

2.3. The evolution of the magnetic potential. Since the magnetic field is divergence- 
free, it can be written as the curl of a magnetic vector potential: 



V x A. 



(2.18) 



The essential feature of constrained transport methods is to update the magnetic vector po- 
tential in some manner, and then by computing from it some discrete form of the curl, to 
obtain a discrete divergence-free magnetic field. Using the relation 

V • (uB - Bu) = V x (B x u) , 
the last row of (2.1) can be written in the form 

B, f |Vx(Bx u) = 0. 
Since B is divergence free, we set B = V x A and rewrite (2.19) as 

V x {A, 4 + (V x A) x u} = 0, 



A t + (V x A) x u= -Vtp, 



(2.19) 

(2.20) 
(2.21) 



where ip is an arbitrary scalar function. Different choices of ip represent different gauge con- 
dition choices, see [26]. 

Our approach is based on the Weyl gauge, which means we take ip = 0. This results in 
an evolution equation for the vector potential is given by 



A, t + (V x A) x u = 0, 
which can also be rewritten in quasilinear form as 

A, t + m (u) A iX + N 2 (u) A „ + N 3 (u) A,. 
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(2.22) 

(2.23) 

(2.24) 



The the flux Jacobian in some direction n = (n 1 ,n 2 ,n 3 ) 



3\T 



M(n, u) := n l Ni + n 2 N 2 + n 3 N 3 



n 2 u 2 + n 3 u 3 



n l u x + n 3 u 3 



n l u x + n 2 u 2 



(2.25) 



has real eigenvalues for all ||n|| = 1, but for any nonzero velocity vector, u, there exist 
directions n, for which M(n, u) does not have a complete set of right eigenvectors. Therefore, 
system (2.23)-(2.24) is only weakly hyperbolic [26]. In order to see this weak hyperbolicity, we 
write the eigenvalues of matrix (2.25): 



A = {0,n • u,n • u}, 



and the matrix of right eigenvectors: 
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Assuming that ||u|| =/= and ||n|| = 1, the determinant of matrix R can be written as 

det(.R) = — ||u|| 3 cos(a) sin 2 (a), 



(2.26) 



(2.27) 



(2.28) 



where a is the angle between the vectors n and u. The difficulty is that for any non-zero 
velocity vector, u, one can always find four directions, a = 0, 7r/2, 7r, and 37r/2, such that 
det(i?) = 0. In other words, for every ||u|| ^ there exists four degenerate directions in which 
the eigenvectors are incomplete. 

3. Outline of the constrained transport algorithm and temporal discretization. 

In our constrained transport method for the MHD equations, wc separate the discretization 
in space and time. For the temporal discretization we use an SSP-RK method of order up to 
three, see e.g. [25, 24]. 

Consider a system of ordinary differential equations of the form 

Q'(t) = C(Q(t)). (3.1) 

One time step of the third order accurate SSP-RK method can be written in the form 

Q« =Q n + AtC(Q n ), 

1 m, 

(3.2) 



Q^ = \q n + \q^ + 1 -mc{qV) 



xn+l 



1 



;Q r 



;Q 



(2) 



-AtC(Q 



(2)] 



In a method of lines approach we obtain systems of ordinary differential equations of the 
form (3.1) after discretizing the spatial derivatives. Here the components of the vector Q(t) 
represent cell average values of the physical quantities in the different grid cells. 

The scmi-discrctc form of (2.1) now has the form 



,(t)=£l(QMHD(t)), 



(3.3) 



where Q M hd (t) represents the grid function at time t consisting of all cell averaged values of the 
conserved quantities form the MHD equations: (p, pu, £,B). The precise form of the spatial 
discretization represented by Ci(Q MHD (t) will be discussed in Section 4. The semi-discrete 
form of the evolution equation for the potential (2.22) has the general form 



?a(*)=£2(<2aW,<2mhd(*)), 
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(3.4) 



where Q a (t) is the grid function at time t consisting of the cell averaged values of the magnetic 
vector potential, A. Note that the evolution of the potential depends on the velocity field, 
which we take to be as given from the solution step of the MHD equation. This is reflected 
in the notation used in (3.4). The spatial discretization of the potential equation will be 
discussed in Section 5. 

A single stage the constrained transport algorithm has the form: 

0. Start with <5™ IHD and Q A (the solution from the previous time step) and Qmhd and 
Q A (the solution from the previous time stage). 

1. Update without regard to the divergence- free constraint on the magnetic field: 

^t!MHD = a VmHD + I t — Of ] Gtf MHD + p AlZ_l I VMHD I I 

Q { ? = « (fe) Q A + (l - « (fe) ) Qt 1] +/3 (fc) At£ 2 (Qj-^Qfc 1 ') , 

where QLhd = {p^,pv^ k \ £^ k \ B^*)) and B( fe *) denotes the predicted value of the 
magnetic field in the first Rungc-Kutta stage. 

2. The magnetic field components of Qm H * d are then corrected by V x Q A . 

B( fc )=VxQ« =► 0L fc H U(p (fe) ,puW,£W,B( fc )). 

For smooth solutions, this procedure with third-order accurate SSP-RK gives an update of 
the grid function QJJhd that is 3 rd order accurate in time. This fact is confirmed numerically 
via test computations done in Section 7. 

Note that in each stage we take one Euler step on the magnetic vector potential equation 
(2.22). In each of these steps we evaluate the spatial discretization operator, £2, at the current 
values of the potential, Qa, and the current velocity values that can be obtained from Q M hd- 
Just as in the constrained transport approaches of [26, 39], during the Euler step on (2.22) we 
view the velocity u as a given function, and view (2.22) as a closed equation for the magnetic 
potential A. 

4. Spatial discretization of the MHD equations. Our spatial discretization of the 
MHD equations is similar to the one used in SHARPCLAW by Ketchcson et al. [27]. However, 
in our method we use a multidimensional polynomial reconstruction. This gives the full order 
of convergence (here up to third order) for smooth nonlinear problems on Cartesian grids and 
reduces grid effects for non-smoothly varying mapped grids. We describe the one-dimensional, 
multidimensional Cartesian, and multidimensional mapped grid spatial discretizations in the 
subsequent three subsections. 

4.1. One-dimensional spatial discretization. We first briefly discuss the one-dimen- 
sional case, where we consider a hyperbolic conservation law of the form 

q.t + /(«),* = 0, (4.1) 

where q : M + x M — > R m is a vector of conserved quantities and / : M. m —> R m is the flux 
function. For smooth solutions we can transform (4.1) into the equivalent quasilinear form 

Q,t + A(q) q iX = 0, (4.2) 

where A(q) := f. q (q) is the flux Jacobian matrix (A : M. m —> R mxm ). We assume that some 
initial data q(0, x) — qo{x) with qo : R — > R m and some appropriate boundary conditions are 
given. 



The cell average of the conserved quantity in grid cell i at time t is denoted by Qi(t), i.e. 

Qi(*)«-^- P +h q(t,x)dx. (4.3) 

Ax Jx. 1 

z 2 

From the given cell average values we compute piecewise polynomial approximations (with 
WENO limiting) of the conserved quantities, using polynomials of degree at most p. For 
smooth solutions these polynomial agree up to order p with the exact solution: 

q(t n ,x):=qi(t n ,x) for x € (x^i ,x i+ i J , 
q i (t n ,x) = q(t n ,x) + 0(Ax p+1 ). 
The reconstructed values from grid cell i at the left and right grid cell interface are denoted 

by 

q + i := lim q~; \X; i +e| and q~ 1 := lim q~i [x, ■ , i — e) . (4.5) 

In semi-discrete form, the method that is used to evolve the cell averages of the conserved 
quantities is given by (see [27]) 



1 

Ax 



Q'M = -— (^~ A <Ii+h + - A+ Aq l _, +AAqA . (4.6) 



The fluctuations A ± Aq i+ i can be computed as in the standard wave propagation algorithm 
of LeVeque [28], or by using the /-wave approach of Bale et al. [1], with the only difference 
being that the left and right values of the conserved quantities used in the wave decomposition 
or /-wave decomposition are obtained from the reconstructed interface values instead of the 
cell average values. The additional term 

f x i+i, r^+i 

AAq t = / A(q~i) q iiX dx « / A(q) q^ x dx, 

J X . 1 J X . 1 

1 2 % 1 

can be computed via Gaussian quadrature of the appropriate order. For the approximation 
of a hyperbolic problem in the flux difference form (4.1), the integral is equal to the flux 
difference: 

2 A{q t ) q iiX dx = f [q~ + ij - / [it-y ■ 

1 
2 

Furthermore, the fluctuations satisfy 

^ + ^-i=/fc)-/t|(C|^)> 

A-Aq l+ i = fl +1 (q-^q+A - f (qT, 



Here /* denotes the interface flux, which is computed by solving a Riemann problem for the 
reconstructed left and right interface values. Thus for hyperbolic problems in divergence form, 
the method can be formulated as a finite volume method of the form 

The Riemann problem can be solved using a variety of solvers; as a matter of practice, in this 
work we use the Roe- type Riemann solvers of the form described in Bale et al. [1]. 
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4.2. Multidimensional Cartesian grids. For simplicity, we consider a two-dimensional 
hyperbolic equation of the form 



or of the quasilinear form 



q,t + f(q), x +g(q), v =0, (4.8) 



q,t + A{q) q tX + B(q) q iV = 0. (4.9) 



Here q : R + x R 2 — >• R m is a vector of conserved quantities, /, g : M. m — > M. m are numerical 
flux functions, and A(q) = f'(q), B(q) = g'{q) are the flux Jacobian matrices. 

In order to get a high-order accurate approximation in space with the method of lines 
approach, we need to reconstruct values on cell interfaces at nodes of a quadrature formula 
that is used to approximate interface fluxes (or fluctuations). In order to do this, we compute 
a piecewise polynomial reconstruction q of the conserved quantities in each grid cell. The 
grid cell average of this reconstructed polynomial agrees with the original cell average (this 
is needed to guarantee numerical conservation). Furthermore, the reconstruction is based on 
a (limited) least squares approach that includes all neighbors that share an edge or a corner 
with the considered grid cell, see [5, 44, 46]. 

The numerical method can be written as a finite volume scheme in the semi-discrete form 
or in the form used by the wave propagation algorithm (compare with [27]) 



At) ~ (U + A 9i i 1 + A~Aq i+i , + AAq 



Ax 
1 

Ay 



B + Aq tl _i+B-Aq^ +k +BAq l3 



(4.11) 



Now Qij(t) is an approximation of the cell average of the conserved quantity at time t in 
grid cell (i,j). Using the piecewise polynomial reconstruction of the conserved quantity q in 
each grid cell, we can evaluate left and right values of the conserved quantity at Gaussian 
quadrature points along grid cell interfaces. These are used to obtain a high order accurate 
representation of the flux: 



(4.12) 



y v o-i y fc=i 

hj, ' +2g *( q ( X > y ^)) dx ™ h^ Ck9 * ( q >^±h<i±h) =:9 W (413) 

J *i-\ fc=i 

The coefficients c^ are the quadrature weights and Xi k and yj k are the quadrature points. In 
an analogous way we can define the fluctuations A ± Aq and B ± Aq along a grid cell interface. 

In order to compute the flux values /* and g* at the Gaussian quadrature points on the 
cell interface, we use a Roe-type Riemann solver of the form described in Bale et el. [1]. An 
advantage of this Riemann solver over the classical Roc solver is that Roe averages can be 
replaced with simpler averages; this feature is very helpful in the case of the MHD equations. 
In any Roe-type solver, we are required to compute right and left eigenvectors of the flux 
Jacobian; for the MHD system one has to be careful about the eigenvector scalings in order to 
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avoid singularities in the eigenvector basis [38] . For the MHD system we use the eigenvectors 
proposed by Barth [6, pages 212-214], which are based on the entropy variables and have near 
optimal scaling properties for physically admissible solution values. 

We have shown here only the 2D spatial discretization. We omit the details of the 3D 
spatial discretization, since the extension to 3D is straightforward. In 3D we use a three- 
dimensional limited piecewise polynomial reconstruction of the conserved quantities and the 
fluxes at grid cell interfaces are computed by integration over a rectangular area (i.e., a surface 
integral instead of a line integral as in 2D). 

4.3. Mapped grids. The spatial discretization can also be extended to two-dimensional 
logically rectangular mapped grids and three-dimensional hexahedral grids. The vertices of 
each mapped grid cell in physical space, (x,y,z), are obtained by mapping the vertices from 
a Cartesian mesh in computational space , (£, r/, £) (sec [11]). 

In three space dimensions, each grid cell can be represented by a trilinear map (also called 
a ruled cell) [10, 47] of the form 

X (C, V, C) = cooo + ciooC + c 010 T] + cooiC + Cno£r? + cioiCC + con??C + cm^(, 

where X = (x,y, z), < £, 77, £ < 1, and Cooo> • • • |Cm € K 3 are vector coefficients. The 
coefficients are determined from the vertices of the hexahedral cell. This representation of 
the grid cell can be used to compute the volume of each grid cell, as well as face areas of all 
the faces, and the unit normal vectors at each face. The two-dimensional case is obtained by 
setting £ — 0. 

Consider a multidimensional hyperbolic system of the general form 

<Z, t + V-F = 0, 

where the columns of the matrix F = [/|g|/i] are the flux functions in the x, y and z direction. 
We consider the semi-discrete form of the finite volume method for a three-dimensional grid 
cell that we denote Cijk- The volume of the grid cell in physical space is denoted by |CV,fc|. 
This volume can be computed using the Jacobian determinant of the trilinear map 



\Cijk\ — 

/o Jo 



ax 



d^drjdC- (4.14) 



Formulas to evaluate this expression are given in [47]. 

The change of the cell average of the conserved quantity in grid cell (i,j, k) is described 
by 

«*<" - w fIL v FdV ■ w £„, F • ""• <" 5 > 

where v is the outward pointing normal vector along the outer boundary of the grid cell. 

At each face of a three-dimensional grid cell, one of the variables £, 77, C is either zero or 
one and the face is represented by a ruled surface. Consider for example the grid cell face that 
corresponds to C, — in computational space (this is the interface with the index (i,j, k — |)). 
In physical space, this interface is described by the bilinear map 

x (£, V) = c oo + cio£ + c Q1 ri + cn£?7, 
11 



which maps a square in computational space to a ruled surface embedded in M 3 . In order to 
further describe this grid cell interface we define vectors 

<9X <9X 

*(1) = -QfT = C l" + C" 7 ?' *(2) = "^ = C m + C "?' 

which are tangent vectors to coordinate lines. The surface metric tensor (aij)i,j=i,2 is defined 
as 

a.j = t ( i) -t U) , i,j,= 1,2, (4.16) 

and a = ana,22 — a,\i(i2\ denotes the determinant of the metric tensor. The unit normal vector 
to the grid cell interface can be computed using 

P(l) x t (2)|| 

Note that with this definition, n$_i „• ^ is an outward pointing normal vector for grid cell 
(i — l,j, k) and an inward pointing normal vector for cell (i, j, fc), compare with the signs of 
the terms in equation (4.19). 

The area element dA on the grid cell interface transforms according to 

dA = y/ad^drj = ||t (1) x t (2) || d^drj. (4.18) 

Note that the surface normal vector n and the determinant a arc functions of £ and r\. Anal- 
ogously we can express the area element and a normal vector at all other grid cell interfaces. 
Using this, we can express (4.15) in computational space by 

Q'iikit) = T^~ 



o Jo 
l r i 



F-n(T,,OVa(v,C)) , -(F-n( V ,CWa( V ,0) 1 ) d v d( 
+ f f {(F-n(e,Ov^O).. 1 -(F■n(tC)^/aJ^0) i ) d^d( 

JO JO I 'ij + ^k x /l 3 2 fe J 

f \(F-n(Z,T,)y/a(Z,T,)) -U-n&vWa&V)) 1 )^dr ] 

(4.19) 

We integrate over grid cell interfaces in computational space using two-dimensional quadrature 
rules. The flux computation normal to the interface is again based on an approximative 
Ricmann solver using the eigenvector decomposition of Barth. 

5. Spatial discretization of the non-conservative magnetic vector potential 
equation. In this section we discuss the spatial discretization that is used to approximate the 
evolution equation for the magnetic potential. For the 2D ideal MHD equations the relevant 
scalar evolution equation for the magnetic potential is (1.1), while in 3D the relevant weakly 
hyperbolic system is (2.23)-(2.24). Just as in the previous section, we introduce the basic 
features of the numerical scheme for a onc-dimcnsional problem. 

5.1. One-dimensional weakly hyperbolic systems. We consider an equation of the 
general form 

q tt + A(x) q, x = 0, (5.1) 
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with g:R + xR^R m and A(x) € R mxm . We assume that system (5.1) is weakly hyperbolic, 
which means that matrix A(x) has real eigenvalues but not always a complete set of right 
eigenvectors. As in the classical hyperbolic case, we wish to construct a method of the form 
(4.6). However, due to the weak hyperbolicity, the fluctuations A ± Aq cannot be approximated 
via a wave decomposition as in [1, 27, 28]. 

Let q be a piecewise polynomial reconstruction of the function q as described in (4.4). 
Furthermore, let A{x) denote a piecewise polynomial approximation of the matrix A(x): 

A(x)=Ai(x) for xe (xi_i,x i+ ij , 
Ai(x) = A(x) + O (Ax p+1 ) . 

The values q x are the reconstructed values at the grid cell interface as defined in (4.5). 

* 3 

Analogously we define 

A7 , := lim A_i \x,- i — e) and A + , := lim A- \x,- i + e) . (5.3) 

Using these reconstructions, we derive a method of the form 



Q'M = -ti M~ A «i+i + ^ +A ?M + -4 A ft 



1 

Ax 
where Qi (t) is the cell average of the quantity q and 



AAqi w lim / A(x) g x dx, (5.4) 

1 2 
/■^»-l/2+£ 

-4 + A&_ 1/2 w lim / ^(x^dx, (5.5) 



1/2 *-<>,, 

r 2: i+i/2 
A~ Aq i+1/2 & lim A{x)q x dx. (5.6) 

e ^° Jx i+1/2 -e 



The integral in (5.4) is approximated by replacing A(x) and q by the approximating polyno- 
mials Ai and gj, i.e. 



AAq t := Ai(x)q~i iX dx. 

J X . 1 



This integral can then be computed using Gaussian quadrature of the appropriate order. 

In the weakly hyperbolic case, the fluctuations A Aq cannot be computed with the wave 
propagation algorithm, since the wave propagation method requires a complete set of eigen- 
vectors. In order to compute these fluctuations, we introduce a regularization of q at each 
grid cell interface of the form 



Ci(^)= *< 



/ 2 \ 




x < x, i — e, 

' 2 


f a; — x i +e \ 




/ 


v - 2 'V 


: x e 


(Si-i-e>3i-i 


CjW 




x > x, i + e, 

' 2 
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where \IJj_i is a path connecting q_i(t) and <z^_i(i). Using this regularization, we first 

2 t 2 * 2 

construct an expression for the sum of the left and right-going fluctuations at the grid cell 
interface x i _i, i.e., we wish to express 

/■X._i+e 

A~Aq t _i +A + Aq i _ 1 « lim / 2 A (a) g x dx. 



This integral is approximated by replacing q with the regularized function q e and by replacing 
the matrix A with the piccewise polynomial approximation A[x), i.e., we set 



A~A qi _i_ + A + Aq^ 



lim / 2 i(x) (a £ _i(i,x)) dx. 



(5.7) 



x— x ._ i +e 



After the substitution s = 

the integral in (5.7) can be written as 

fX ._ i +e 



and by defining 0(s) := .Tj_i — e + s{x i _ i +e— (xj_ i — e)), 

(5.8) 



A{x)[ql_,{t,x)\ dx = A(<f,{a)) [%_ h (a,t)) ds 

x. i-e V 2 y >x j 



In order to resolve this integral, we choose the simple straight-line path 

*»-§ = ^~ i + s (gti - «~_i ) i 
and obtain 






l(^( S ))(^_i( S ,t)) ds= f A Ms)) ds (q+_ ,{t)-q-_ ,{t) 

IQ V * J ,8 Jo V 2 2 

The integral on the right-hand side can be split into two pieces: 

1 .1/2 f l 

A{4>(s))ds= A i _ 1 (<p(s))ds+ A\{<f>{8))ds, 

Jo J 1/2 



(5.9) 



(5.10) 



(5.11) 



which corresponds to a division into pieces on the left and right of the discontinuity x i _ i . We 
then take the limit e — > and get: 



lim 



1/2 



1 



1 



A i _ 1 (^( S ))d S + / Ai($(s))ds = -A- k + -4+ =: A 



1/2 



2 l - 



Thus we obtain 



A A 9i _i+^ + A^_i =A\ 9 Jg+i-^i). 



(5.12) 



(5.13) 



Following Castro et al. [13] and using relationship (5.13), we are now able to define the 
left and right-going fluctuations. With partial knowledge of the eigenstructure of A\^ (i.e., if 
we know the eigenvalues or an estimate of the largest absolute value of the eigenvalues) we 
can define 



A Aq i+1/2 = 



1 



A 



*« + l/2 



Oii+1/2^ 



q i+h- q i + i 



(5.14) 



generalized Rusanov flux 
14 



and 



A + Aq t _ l/2 



1 



A 



, +ai_i/ 2 ] 

i-1/2 ' 



(5.15) 



generalized Rusanov flux 

where I £ ]]j mxm j s the identity matrix. Here a is a positive number with 

|A fc | < a, for k — 1, . . . , m, 
where A fc represents the fc th eigenvalue of Al , see [13]. 

Another possibility is to define the fluctuations without using any knowledge of the eigen- 
structurc of A L : 



A Aq i+1/2 



1 



2A\ 



Ax 
'At 



^1 II 
A.: i A| * 



*H- 



H+i 



and 



A + Aq^ 1/2 = 



2A 



*!-* At 



A *i+^UI 



Aa 



9, 



(5.16) 



(5.17) 



These fluctuations are derived from the generalized FORCE scheme [12, 13], which is a convex 
combination of the generalized Lax-Friedrichs and the generalized Lax-Wendroff scheme. 

5.2. Multidimensional Cartesian grids. We now consider the approximation of a 
weakly hyperbolic system of the form 



with q : IS 
the form 



q,t + A(x, y) q iX + B(x, y) q. y = 0, (5.18) 

l , A(x,y), B(x,y) £ u mxm . i n the semi-discrete case, the method has 



Q'n(t)=- 



Ax 
1 

Ay 



A + Aq t 



B + Aq 



*,J- 



:j + A Aq i+iJ +AA gij 
h +B-Aq ij+i +BA qij 



(5.19) 



As in the case of hyperbolic systems (see Section 4.2), we compute a multidimensional piece- 
wise polynomial reconstruction q of q using a least squares approach. Furthermore, A(x, y) 
and B(x, y) are piecewise polynomial approximations of the matrix valued functions. Using 
these reconstructions we define 



AAq i:j := 
BA qij := 



1 
A^ 

1 
Ax 



Vi+l f x i+\ - 

/ A ij (x,y)q ijjX dxdy, 

Vi_l •'z. i 

3 2 2 

\ B lj (x,y)q i j tV dxdy. 



(5.20) 



The fluctuations A and B arc defined in analogy to (5.14) and (5.15) with the only difference 
that we integrate over a grid cell interface and then compute the average value. For example, 
we compute 



A~Aq. 



i+%j 



i_ r v i+i i 



A\ T (y) - a,, i(v)I 
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q t + y(y)-Q- +lJ (y))dy. (5.21) 



We evaluate this integral using Gaussian quadrature formulas. This means that we need to 
compute the matrix ^4 at the nodes of the quadrature formula. The other fluctuations 

are computed in an analogous way. 

5.3. Mapped grids. We now discuss the discretization of a weakly hyperbolic system 
of the form (5.18) on a logically rectangular mapped grid. For simplicity we restrict our 
considerations to the two-dimensional case. The extension to the three-dimensional case is a 
combination of the concepts discussed in this section and those of Section 4.3. 

Let Cij denote a grid cell in physical space and |CV,| the area of the grid cell. Assume that 
the grid cell in physical space is obtained by a mapping of the unit square in computational 
space. The change of the cell average of the quantity q in grid cell Cij is described by the 
ODE: 

Qij{t) = -TTTlUi-ij A + Aq t _i j + A tj q ijjX + £ l+ i 3 A~Aq i+ i i 

|L/ ^ |V " " " (5.22) 

+ t {j _ i B+Aq itj _ i + B l3 q ijiV + t ij+ i B-Aq iJ+l " 

where £ i+ i and ^,- + i denote the length of the grid cell interfaces in physical space. The 
integrals over the grid cell arc defined as 

Aj qij, x = // A-ij (x, y) qi 3 .x (x, y) dx dy 



\/aij(£, V) Aij(x(£, r)),y(€, V)) QijA x (^ V),y(^i V)) dS, dr/, (5.23) 



Jo 
1 r l 



o Jo 



where we use piecewise polynomial reconstructed functions q and a piecewise polynomial 
representations of the matrix valued functions. For a planar two-dimensional mapped grid 
cell, the area element \fa is defined by 

Va = | {c\ + c{ lV )(c 2 Q1 + c?^) - (cji + ciiOC^o + cur?) 

The coefficients Cio, • • • , Cu £ R 2 are special cases of those detailed in §4.3, and can be obtained 
from those coefficients by setting £ and the third component of X equal to zero. 

In order to define the fluctuations A ± Aq and B ± Aq, we introduce n i± i and r^ ± i , the 
normal vectors at the interfaces (i± ^ j) and (ij ; ± g)- We define A L by replacing the matrix 
A in the one-dimensional formula of Section 5.1 by A = n 1 A + n 2 B. Now we integrate the 
fluctuations over each grid cell interface using again the piecewise polynomial reconstructed 
values. For example we compute 



*'**»*>■- J 2 



M 



A\y {v) - ",+iW 1 

» -t- — -i z 



<lt +hj {ri)-q- +lj {r } ))dr t . (5.24) 



All the other fluctuations are computed analogously. 



5.4. Limiting with respect to the derivative. As first pointed out by Rossmanith 
in [39], a special limiting is needed for the update of the magnetic potential in order to 
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avoid unphysical oscillation in derivatives of the potential and thus in the magnetic field. He 
proposed an extension of the limiting mechanism used in the wave propagation algorithm. 

Here we describe a limiting strategy that can be used in combination with the update 
of the potential that is proposed in this paper. The strategy of limiting with respect to the 
derivative will be described for the ID advection equation 

q,t+ m.x=0, (5.25) 

where q : R x M + — > R is continous but not necessarily continously diffcrcntiable. The 
advection speed u G M. is assumed to be constant and we restrict to the case u > 0. We use 
the method introduced in Section 5.1 together with the third order accurate SSP-RK method 
to update the advcctcd quantity q. 

In order to limit the solution we add numerical viscosity 1 to the problem, i.e. instead of 
(5.25) we approximate an advection diffusion equation of the form 

q,t + uq,x = e(ar) q, x ,x- (5.26) 

This scalar advcction-diffusion equation can also be written as a system of first order equations 

q,t + uq yX = e(x) d, x , (5-27) 

d-q, x = 0. (5.28) 

Usually e(x) is chosen such that a thin layer is added near discontinuities in the solution itself 
while keeping the high order reconstruction away from the discontinuity (e.g., see [34]). We 
will change the definition of e in such a way that additional viscosity is added near jumps in 
the derivative of the solution. 

The artificial viscosity must satisfy the following requirements: 

1. It should be small enough to satisfy the stability constraint 

e{x)At 1 
(Ax) 2 " 2 

of our explicit time-stepping method. 

2. It should be large enough to avoid spurious oscillations in the derivative due to Gibbs 
phenomenon. 

3. The artificial viscosity should not degrade the acurracy of the scheme for smooth 
solution structures. 

To be able to handle all these requirements we define the viscosity in the following way: 

e(x) = r]a, (5.29) 



with 



71 = ' 2 ^a7-' (5 ' 30) 



and 

_ Ji[l + sin(7rAS-f)] S > a u , 



(5.31) 
S < an, 



1 In the context of the evolution of the magnetic potential in ideal MHD, the term numerical viscosity 
really refers to numerical resistivity. 
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where 



A 



&ik 



i k 



((Ax)* + z ik y 



S = max((7ii_i,o-ii+i), AS=\S-<ru\, 
S 4fe = (q }XtX (k)(Ax) 2 ) 2 \ k e (* - 1,*,* + 1). 



(5.32) 



The parameters are set to Aj 



A 



»t+l 



1, A, 



1000 and 



4 in the test below. Note 



that we use a double index on the parameters A and a; the logic behind this notation is that 
the first index denotes the current cell that is being considered and the second index denotes 
the cell to which the current cell is being compared. r\ is the maximum amount of viscosity 
that could be added and a is a smoothness indicator, with values between and 1, such 
that a approaches 1 on very steep gradients in the derivative (i.e., steep second derivative), 
and a is in smooth regions of the derivative. The idea behind this smoothness indicator 
is that we do not take any viscosity (large linear weight An), unless the second derivative in 
one of the neighbouring cells is much higher than in cell i, in which case we will add some 
viscosity depending on the difference of the smoothness measure an and the maximum of the 
smoothness measures in the neighouring cells (an-i,an+i). This means that the parameter 
Aj; controls how much the second derivatives have to differ, before adding any viscosity. We 
found that An G [10 1 , 10 3 ] seem to give good results. From the form of the viscosity, it is clear 
that we do not obtain a further time step restriction due to the viscosity limiting. Taking the 
time step restriction into account, one could also choose r\ = O(Ax) empirically. 

We now consider the ID advection example that was first introduced in [39]. Consider 
(5.25) with u = 1 on the interval < x < 1 and double periodic boundary conditions. The 
initial data is the following piecewise linear function: 



9(0, x) 





(x - 0.25)/0.075 

2 

(0.75 - .t)/0.075 





if x < 0.25, 

if 0.25 < x < 0.4, 

if 0.4 < x < 0.6, 

if 0.6 < x < 0.75, 

if 0.75 <x. 



Note that q(0, x) is continuous, but its first derivative is discontinuous. The challenge is to 
control oscillations in both q{t,x) and q tX (t,x). 

We compare the solution and its derivative as computed by two different limitcrs: (a) a 
standard WENO limitcr that is essentially non-oscillatory in q(t,x), but not necessarily in 
q, x (t,x), and (b) the proposed limiter as described above. For further details on the WENO 
limiting see Canestrclli ct al. [12], Titarcv ct al. [44] and Tsoutsanis et al. [46]. Third order 
accurate average values of the derivatives q tX (t,x) are computed using 



1 

Ax 



q, x (t,x)dx 



1 
Ax 

1 
Ax 



q(t,Xi 



Q\t,Xi_i 



^W + ^-^CjW + ^W) 



(5.33) 



We use a CFL number of 0.7 in both simulations. The solution and its derivative at time 
t = 1 (i.e., after one revolution) are shown in Figure 5.1, where Panel (a) is the WENO limitcr 
and Panel (b) is the proposed limiter. The approximation of q(l,x) looks very similar in 
both cases. However, there are obvious differences in the computed results for q iX (t, x); the 
proposed limiter avoids spurious oscillations in q yX (t, x). In comparison with the TVD scheme 
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Fig. 5.1. The solution to the ID advection equation and its derivative on a periodic domain after one 
revolution using the 3 rd order algorithm. Panel (a) shows how standard WENO limiting performs on the 
solution, and Panel (b) shows how the limiter proposed in this work performs. 



of [39] the limiting introduced here (which can also be used in a one-step scheme) seems to 
be slightly more diffusive. This however is a general problem if a method of lines approach is 
used instead of a Lax-Wendroff type method. 

In the multidimensional case, we compute the smoothness indicator S by taking all neigh- 
boring grid cells into account that share a face, an edge or a corner point with the considered 
grid cell. For the computation of £, the second derivative of q is replaced by the Laplacian and 
(Ax) 2 is replaced by (|Cjj|) in the two-dimensional case, and (|CV,-fc|) in the three-dimensional 
case. As characteristic length in the computation of r\ and a (Ax in ID), we use the length 
of the smallest edge of the considered grid cell in the multidimensional case. Denoting the 
characteristic length by As, we use 77 = 0.5As in all MHD computations below and Aj j_i = 1, 
Ajj+i = 1, and A« = 10 3 (see formulas in (5.32)). 



6. The discretization of V x A. During each stage of the constrained transport al- 
gorithm described in Section 3, we obtain the corrected magnetic field from taking the curl 
of the magnetic potential. In our previous work on CT schemes [26, 39], the cell average 
values of the magnetic field were computed using centered finite difference approximations. 
In order for the scheme developed in this work to be third order accurate both on Cartesian 
and mapped grids, we must generalize the central finite difference approach. We describe this 
generalization in this section. 
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We wish to compute a discrete version of the curl of the magnetic vector potential: 

B = VxA=[A 3 y -A% A] z -A% A% ~ A] y ] T . (6.1) 

In particular, we wish to compute a third order accurate estimate of the cell average of B. For 
example, using the divergence theorem, the cell average of the first component of the magnetic 
field, B , in grid cell (i,j, k) can be expressed in the form 



i 




(AV - A 2 v 3 ) dA, (6.2) 

ldC ijk 

where v — (y\ v , u ) G M 3 is the outward pointing unit normal vector along the boundary 
dCijk of the considered grid cell. The integrals over all the faces of a grid cell can be evaluated 
in comptational space as discussed in Section 4.3. Using the definition of the surface normal 
vector n and the determinant a of the mapping as defined for the grid cell interfaces in Section 
4.3, we obtain 

Bl jk = r^—, E { ± / / (V^O - A 2 n 3 { v X)) V^O) ^ t <MC 

\ C ijk\ ±t£,- L JO JO V >i±\,3M 

iffteo-AVlcc))^ x dCdC ( 6 - 3 ) 

JO JO V /i,j±±,k 

±f f ((A'V(e,r / )-A 2 n 3 (^,r ? )) v ^(^) i d£dri\. 

Jo Jo v /*,j,fe±2 > 

In order to evaluate the surface integral in (6.3), we first replace the integrals with a 
Gaussian quadrature rule of the appropriate degree of precision, then we reconstruct A 2 and 
A 3 on the grid cell boundaries at each of the quadrature points using the piecewise polynomial 
reconstruction of A. To be more precise, at each of these quadrature points we take the value 
of A to be the average of the reconstructed values from the two grid cells that share a face. We 
note that this discretization leads to a conservative update of the magnetic field. For a smooth 
magnetic potential the average of the reconstructed values agrees with the correct interface 
value of the magnetic potential up to the order used in the reconstruction. The computation 
of the cell averages of B 2 and B 3 can be done in an analogous way using 

B U = T7^$ i Al " 3 A ^) dA 
1 



\Cijk\ JJdC iik 



ijk 



\Cijk\ jjdc, 

Claim 6.1. The constrained transport method as described in this work locally conserves 
(and therefore also globally conserves) the magnetic field, B, from one Runge-Kutta stage to 
the next. 
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Proof. Without loss of generality consider a single Euler step (i.e., a single stage in the 
SSP time-stepping scheme) on the magnetic field: 



B n+1 = V x A" +1 = V x (A™ + A*E n ) 

= B" + At V x E" = B™ + At V • ( [e ijk E k ] 



,-M (6-4) 



The divergence of the tensor [eijkE k ] is computed via surface integrals of the form (6.3), 
where A is replaced with E. The values of E = (V x A) x u arc computed from constrained 
transport as described in Section 5 and then averaged onto appropriate Gaussian quadrature 
points on each face. In particular, note that the same values of E are used to update the 
cell average of V • ([e^fc-E*]) on either side of the face. Therefore, we are guaranteed that the 
discrete magnetic field is locally conserved over each Euler time step. D 

7. Numerical Results. Several numerical examples on Cartesian and mapped grids 
are shown in this section. These examples are used to both verify the third-order accuracy in 
space and time, as well as the shock-capturing ability of the scheme proposed in this work. 

7.1. Smooth Alfven wave problem. We consider two variants of the so-called smooth 
Alfven wave problem: a (1) 2.5D variant and (2) 3D variant. Note that the term 2.5D here 
refers to the case of a two-dimensional computational domain, but with vector unknowns u 
and B that have three non-trivial components. In 2.5D for the magnetic potential we solve 
equation (2.23) with A z — using the method described in Section 5. 

7.1.1. 2.5D problem. In this case the analytical solution consists of a sinusoidal wave 
propagating at constant speed in direction n = (cos 0, sin 0, 0) without changing shape. We 
use 4> = arctan(0.5) and solve in the domain (x,y) € [0, (cos</>) -1 ] x [0, (sin^) -1 ] with double 
periodic boundary conditions. The initial values of all components are described in [26]. 
Similar problems were also considered by Rossmanith [39] and Toth [42]. 

In Table 7. 1 we show the results of a numerical convergence study of the new constrained 
transport algorithm on a 2D Cartesian grid. We obtain full third order convergence rates in 
all conserved quantities, as well as all the magnetic vector potential components. 

We have also computed this test problem on a mapped grid, which is a scaled version of 
a grid from Colella et al. [15]. The mapping has the form 

*(*,*) -(j|) +M n(^)-=(^)(;), (T.I) 

where (x Cl y c ) are the coordinates in computational space, /3 € M. is a parameter which deter- 
mines the smoothnes of the grid, and L, M describe the length of the domain in the x and y 
direction, respectively. Here we use j3 — 0.1, L = (cos^) -1 and M — (sin^)^ 1 . The mapped 
grid is shown in Panel (a) of Figure 7.1. Table 7.2 confirms the third order convergence rate 
also for the mapped grid computation. Note that the error on the mapped grid is only slightly 
larger than the error on the Cartesian grid. However, we note that the least squares approach 
that we used for the piecewise polynomial reconstruction leads to the correct order only on 
grids that do not have highly stretched cells (e.g., see Petrovskaya [35]). 

7.1.2. 3D problem. Next we have performed a convergence study for the 3D smooth 
Alfven wave problem on a Cartesian grid. The results, which again confirm the third order 
convergence rate, are presented in Table 7.3. The initial values and the computational domain 
for this problem are described in [26]. 
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P 


pu 1 


pu 2 


pu 6 


£ 


32 x 64 


5.128 x 1CT 4 


2.413 x 10~ 4 


6.059 x 10~ 4 


5.142 x 10~ 4 


1.301 x 10~ 4 


64 x 128 


6.495 x 10" 5 


3.097 x 10" 5 


7.548 x 10" 5 


6.469 x 10" 5 


1.653 x 10" 5 


128 x 256 


8.150 x icr b 


3.924 x 10" b 


9.412 x 10" b 


8.099 x 10" b 


2.075 x 10" b 


256 x 512 


1.020 x 10" 6 


4.933 x 10- 7 


1.176 x 10" 6 


1.014 x 10" b 


2.599 x 10~ 7 


EOC 


2.998 


2.992 


3.000 


2.998 


2.997 





B 1 


B 2 


B 3 


32 x 64 


2.877 x 10~ 4 


5.754 x 10~ 4 


5.123 x 10" 4 


64 x 128 


3.583 x 10" 5 


7.166 x 10" 5 


6.437 x 10" 5 


128 x 256 


4.464 x 10~ b 


8.928 x 10" b 


8.057 x 10" b 


256 x 512 


5.581 x 10" 7 


1.116 x 10" 6 


1.008 x 10" 6 


EOC 


3.000 


3.000 


2.999 






A 1 


A 2 


^l 3 


32 x 64 


3.583 x 10- b 


6.921 x 10~ 5 


1.053 x 10- 4 


64 x 128 


4.489 x 10~ b 


8.628 x 10~ b 


1.327 x 10" 5 


128 x 256 


5.626 x 10~ 7 


1.076 x 10~ b 


1.661 x 10" b 


256 x 512 


7.041 x 10~ 8 


1.344 x 10" 7 


2.077 x 10" Y 


EOC 


2.998 


3.000 


2.999 



Table 7.1 

Convergence study of the 2.5D smooth Alfven wave problem on a Cartesian mesh. The tables show the 
Li-error at time t = 1 in the different physical quantities computed using the constrained transport algorithm. 
The experimental order of convergence (EOC) is computed by comparing the error for the two finest grids. 
All simulations are performed using a CFL number of 0.5. 




(a) 



(I)) 






0.6 0.8 



Fig. 7.1. Shown in this figure are the mapped grids used in (a) the convergence study of the smooth 
Alfven wave problem and (b) the cloud-shock interaction problem. 



7.2. 2.5D shock tube problem on mapped grids. Next we consider a ID shock tube 
problem on a scaled versions of the 2D mapped grid given by (7.1). This example demonstrates 
the performance of the new CT method on shocks not moving along grid lines. 

The computational domain is the square [—0.7, 0.7] x [—0.7, 0.7]. In order to avoid inaccu- 
racies at the boundary of the domain, we restrict the mapping to the area —0.6 < x c , y c < 0.6. 
The rest of the mesh is Cartesian. We set L = M = 1.2 and use different values of /?. 
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P 


pui 


pu 2 


pu 3 


£ 


32 x 64 


1.975 x 10~ 3 


7.167 x 10~ 4 


1.223 x 10~ 3 


7.711 x 1CT 4 


5.230 x 10~ 4 


64 x 128 


2.631 x 1CT 4 


9.756 x 10" b 


1.557 x 10" 4 


9.745 x 10~ 5 


6.955 x 10" 5 


128 x 256 


3.333 x 10- b 


1.245 x 10" 5 


1.953 x 1(T 5 


1.222 x 1CT 5 


8.820 x 10" b 


EOC 


2.981 


2.970 


2.996 


2.996 


2.980 





B 1 


B 2 


B 3 


32 x 64 


ami x io~ 4 


9.735 x 10~ 4 


9.483 x 10~ 4 


64 x 128 


6.014 x 10" 5 


1.213 x 10" 4 


1.219 x 10" 4 


128 x 256 


7.511 x 10" 6 


1.513 x 10" 5 


1.534 x 10~ 5 


EOC 


3.001 


3.002 


2.990 






A 1 


A 2 


A* 


32 x 64 


7.355 x 10~ 5 


1.260 x 10~ 4 


1.749 x 10~ 4 


64 x 128 


9.358 x 10" b 


1.595 x 10" 5 


2.216 x 10" 5 


128 x 256 


1.175 x 10" b 


2.000 x 10" b 


2.776 x 10" b 


EOC 


2.994 


2.996 


2.997 



Table 7.2 
Convergence study of the 2.5D smooth Alfven wave problem on a mapped grid. The tables show the 
Li-error at time t = 1 in the different physical quantities computed using the constrained transport algorithm. 
The experimental order of convergence (EOC) is computed by comparing the error for the two finest grids. 
All simulations are performed using a CFL number of 0.5. 





P 


PU\ 


pu 2 


PU-3 


£ 


16 x 32 x 32 


2.055 x 10~ 3 


9.992 x 10~ 4 


1.716 x 10- a 


1.911 x 10- 3 


5.133 x 10~ 4 


32 x 64 x 64 


2.684 x 10" 4 


1.343 x 10" 4 


2.213 x 10" 4 


2.532 x 10" 4 


6.611 x 10" 5 


64 x 128 x 128 


3.419 x 10" 5 


1.679 x 10" 5 


2.764 x 10" 5 


3.110 x 10" 5 


8.701 x 10-° 


EOC 


2.973 


3.000 


3.001 


3.025 


2.926 





L\ Error in B 1 


L\ Error in B 2 


L\ Error in B A 


16 x 32 x 32 


1.231 x 10- 3 


1.805 x 10- a 


1.778 x 10~ 3 


32 x 64 x 64 


1.581 x 10" 4 


2.304 x 10" 4 


2.296 x 10" 4 


64 x 128 x 128 


1.957 x 10" 5 


2.869 x 10" 5 


2.818 x 10" 5 


EOC 


3.014 


3.006 


3.026 





L\ Error in A 1 


L\ Error in A 2 


L\ Error in A 6 


16 x 32 x 32 


1.588 x 10~ 4 


2.933 x 10~ 4 


2.869 x 10" 4 


32 x 64 x 64 


2.340 x 10~ 5 


3.642 x 10- 5 


3.671 x 10- b 


64 x 128 x 128 


2.550 x 10" b 


4.780 x 10" b 


4.617 x 10" b 


EOC 


3.198 


2.930 


2.991 



Table 7.3 

Convergence study of the 3D smooth Alfven wave problem on a Cartesian grid. The tables show the 
Li-error at time t = 1 in the different physical quantities computed using the constrained transport algorithm. 
The experimental order of convergence (EOC) is computed by comparing the error for the two finest grids. 
All simulations are performed using a CFL number of 0.6. 



2:-! 



The Riemann initial data have the form 

(p, u\u 2 , u 3 ,p, B\B 2 , B 3 )(0,x 

_2 

/Sir ' \/3ir ' v / 3ir 



1.08, 1.2, 0.01, 0.5, 0.95, -§=, -#=, -%=) if x < 0, (7.2) 



1, 0, 0, 0, 1, -!=, -4=, -JL ] if x > 0. 

' ' ' ' ' V47T ' V4ir ' V4ir J ~ 

As the initial condition for the magnetic potential we use 

(A 1 , A 2 , A 3 ) (0, x) = (0, xB 3 , yB 1 - xB 2 ) . (7.3) 

We utilize zeroth order extrapolation on all boundaries for the MHD variables and linear 
extrapolation for the magnetic potential. Shown is a comparison of scatter plots for the 
magnetic field components along the x-axis. 

In Figure 7.2, we show results for the magnetic field components using the new CT 
method. Here we compare results for differnt grids, namely a Cartesian grid and two different 
versions of the mapped grid obtained by setting the parameter (3 to either 1/50 or 1/15. 
For all simulations we used grids with 200 x 200 mesh cells. Note that for the mapped grid 
computations the solution structure, moving only in x-direction, is not aligned with the grid. 
It has been observed previously that this leads to unphysical oscillations, see [26, 32]. A 
numerical convergence study for B 1 on mapped grids with different resolution showed that 
the error in the L\ norm decreases with a rate of about 1/3. By introducing more numerical 
viscosity (i.e. decreasing the value of \u) these oscillations can be slightly reduced. However, 
they could not be avoided. Here we show results for \n — 10 3 . The solid lines in these plots 
are obtained by computing solutions of the one-dimensional Riemann problem on a very fine 
equidistant mesh. 

7.3. Cloud-shock interaction problem. Next we present numerical results of the 
proposed method for a cloud-shock interaction problem. Again we consider 3D and 2.5D 
variants of this problem. 

7.3.1. 3D problem. The initial conditions consist of a shock that is located at x = 0.05, 
with 

(p, u\u 2 , u 3 ,p, B\B 2 , B 3 )(0,x) 

(3.86859, 11.2536, 0, 0, 167.345, 0, 2.1826182, -2.1826182) if x < 0.05, (7-4) 
(1, 0, 0, 0, 1, 0, 0.56418958, 0.56418958) if x > 0.05, 

and a spherical cloud of density p = 10 with radius r — 0.15 centered at (0.25,0.5,0.5). 
The cloud is in hydrostatic equilibrium with the fluid to the right of the shock. The initial 
conditions for the magnetic potential are given by 

. . _ f (2.1826182 y, 0, -2.1826182 (x - 0.05)) T if x < 0.05, , . 

^ ' X ^~ | (-0.56418958 y, 0, 0.56418956 (x - 0.05)) T if x > 0.05. 

The computational domain is the unit cube. Inflow boundary conditions are used at the left 
side and outflow boundary conditions are used at all other sides. Without the constrained 
transport step the method would fail to compute the solution structure. In Figure 7.3 we 
present results of a computation with the three-dimensional method proposed in this paper. 
Compare also with [26], where this problem was computed with our previous approach. 
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Fig. 7.2. Shown are the magnetic field components for the shock tube problem at t = 0.2. The left column 
shows the results for a Cartesian grid computation; the middle column shows results on a mapped grid which 
is a small perturbation of the Cartesian grid (using f3 = 1/50J; the right column shows results on a mapped 
grid which is a larger perturbation of the Cartesian grid (using (3 = 1/15,). 



7.3.2. 2.5D problem. We also studied the cloud-shock interaction problem in the 2.5D 
case. The initial values are taken from the 3D test problem by setting z = 0.5. In the 2.5D 
case we show results for B 3 as computed by the new constrained transport algorithm that 
updates all three components of the magnetic field. The results compare well with those 
obtained by our previous approach [26] and also with the results from the 2D unsplit method 
of [39] in which only B 1 and B 2 are updated by a constrained transport method. We compute 
the solution both on a Cartesian grid as well as on a mapped grid. For the mapped grid 
computation we use a grid of the form shown in Panel (b) of Figure 7.1. This is a small 
modification of a grid with circular inclusion discussed in [11]. The initial cloud position is 
inside the circular region of the mapped grid. The results on both grids compare well, although 
the solution structure is slightly sharper resolved on the Cartesian mesh (see Figure 7.4). 



8. Conclusions. In this work we developed a finite volume method for solving the 2D 
and 3D ideal MHD equations on both Cartesian and logically Cartesian mapped grids. The 
scheme is a method of lines that is based on a finite volume discretization in space coupled 
with a strong-stability-preserving Runge-Kutta time stepping method. By using this method 
of lines discretization, we were able to construct a third-order accurate algorithm for the MHD 
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Density at time t=0.06 




1 1 



Fig. 7.3. This figure shows the density at time t = 0.06 on a three-dimensional Cartesian grid using 
150 X 150 X 150 mesh cells. For this simulation we used a 2 nd order version of the proposed algorithm. 



B„ at time 0.06 



B„ at time 0.06 



(a) 




(b) 




Fig. 7.4. The 3 rd component of the magnetic field at time t = 0.06 (a) on a 256 X 256 Cartesian grid 
and (b) on the mapped grid shown in Figure 7. 1 (right) using the 3 rd order version of the algorithm. Here 
the 2.5 dimensional problem was implemented in such a way that all three components of the magnetic field 
were updated in the constrained transport step. 



equations. In order to control errors in the divergence of the magnetic field we developed a 
novel constrained transport approach that couples to the finite volume discretization. Our 
constrained transport methodology allowed us to overcome two important difficulties: (1) the 
evolution equation for the magnetic potential is only weakly hyperbolic, and (2) standard 
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limitcrs applied to the magnetic potential do not adequately control unphysical oscillations in 
the magnetic field. 

The method described in this work is based on the following procedure that is carried out 
during each Runge-Kutta stage: 

1. Update the MHD variables without regard for controlling discrete V • B errors. The 
updated magnetic field in this step is the predicted magnetic field. This step is carried 
out with the finite volume schemes described in Section 4. 

2. Update the magnetic vector potential by solving a weakly hyperbolic vector transport 
equation. This transport equation is simply the induction equation with ideal Ohm's 
law, but written in potential form using the Weyl gauge. This step is carried out with 
the non-conservative finite volume schemes described in Section 5. Special artificial 
resistivity limiters described in Section 5.4 are used in this step to simultaneously 
control unphysical oscillations in the magnetic potential and the magnetic field. 

3. Define the corrected magnetic field as the result of taking the discrete curl of the 
updated magnetic vector potential. The precise form of the discrete curl operation is 
described in Section 6. 

The resulting scheme was applied to several 2.5D and 3D test cases on both Cartesian and 
mapped grids. These test cases demonstrated two important features: (1) we are able to obtain 
full third-order accuracy on smooth problems and (2) we are able to accurately capture shock 
waves. 

Finally we note that although there are some free parameters required in the definition of 
the artificial resistivity added in the magnetic vector potential equation update, we have found 
parameter values that seem to robustly work for a large class of problems. In all the simulations 
presented in this work we have stuck to these same parameter values. Experimentation with 
the Xu parameter in (5.32) shows that one can slightly increase or reduce the amount of 
artificial resistivity, but that for a broad range of A^ the results are qualitatively the same. 
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